      program bernstein_test
!     Solution of two-point BVP for ODE by collocation method with Bernstein polynomials
      use prec
      use bernstein
      implicit none

      real(dp) :: a,b  ! interval
      integer,parameter :: n = 4 ! Order of polyniomial
      real(dp), dimension(0:n) :: nd ! nodes (values in interval [a,b]), there are n+1 nodes!!, where n is order of Bernstein polinomial 

      integer :: i
      real(dp), parameter :: e = EXP(1.0_dp) ! Euler constant
      real(dp) :: tmp,tmp2,t ! sranje za testove

!     Define solutio interval [a,b]
      a = 0.0_dp
      b = 2.0_dp 


!     Form the vector of nodes in [a,b] interval
      write(*,*)
      write(*,'(a)') 'Equidistant collocation points on an interval:'
      nd = bernstein_basis_nodes(n,a,b)
      write(*,*) (nd(i), i=0,n,1)


!    Test Basis fnction and derivatives
     write(*,*)
     write(*,'(a)') 'Writing test results, if correct all results should be     0.0000000000000000  : '
     t=0.5_dp
!    Test I
     tmp = bernstein_basis_fun_eval(2,9,0.0_dp,10.0_dp,t)
     tmp2 = 9.0d0*(10.0d0-t)**7*t**2/250e+6
     write(*,*) tmp-tmp2 ! Trebalo bi da je    0.0000000000000000 !
!    Test II
     tmp = bernstein_basis_fun_eval(1,9,0.0_dp,10.0_dp,t)
     tmp2 = 9.0d0*(10.0d0-t)**8.0d0*t/1e+9
     write(*,*) tmp-tmp2 ! Trebalo bi da je    0.0000000000000000 !  
!    Test III
     tmp = bernstein_basis_fun_derivative(1,2,9,0.0_dp,10.0_dp,t)
     tmp2 = real(9)/10.0_dp*(bernstein_basis_fun_eval(1,8,0.0_dp,10.0_dp,t)-bernstein_basis_fun_eval(2,8,0.0_dp,10.0_dp,t)) ! eq. 5 Bhatti Bracken J.Comp.App.Math. 205 (2007) pp. 272-280
     write(*,*) tmp-tmp2 ! Trebalo bi da je    0.0000000000000000 !
!    TEST IV
     tmp = bernstein_basis_fun_derivative(1,1,9,0.0_dp,10.0_dp,t)
     tmp2 = real(9)/10.0_dp*(bernstein_basis_fun_eval(0,8,0.0_dp,10.0_dp,t)-bernstein_basis_fun_eval(1,8,0.0_dp,10.0_dp,t))
     write(*,*) tmp-tmp2 ! Trebalo bi da je    0.0000000000000000 !


!    Symmetry test - if derivative is even order - these are equal, if der. is odd they are anti-symmetric
     write(*,*)
     write(*,'(a)') 'Symmetry test - values should have same absolute values but opposite signs  : '
     tmp = bernstein_basis_fun_derivative(3,2,9,0.0_dp,1.0_dp,0.8_dp)
     write(*,*) tmp
     tmp2 = bernstein_basis_fun_derivative(3,7,9,0.0_dp,1.0_dp,0.2_dp)
     write(*,*) tmp2

     stop
     end program
